Systematic uncertainties in the determination of the local dark matter density 
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A precise determination of the local dark matter density and an accurate control over the corre- 
sponding uncertainties are of paramount importance for Dark Matter (DM) searches. Using very 
recent high-resolution numerical simulations of a Milky Way like object, we study the systematic 
uncertainties that affect the determination of the local dark matter density based on dynamical 
measurements in the Galaxy. In particular, extracting from the simulation with baryons the orien- 
tation of the Galactic stellar disk with respect to the DM distribution, we study the DM density 
for an observer located at ~8 kpc from the Galactic center on the stellar disk, po. This quantity 
is found to be always larger than the average density in a spherical shell of same radius po, which 
is the quantity inferred from dynamical measurements in the Galaxy, and to vary in the range 
Po/po = 1.01 — 1.41. This suggests that the actual dark matter density in the solar neighbourhood 
is on average 21% larger than the value inferred from most dynamical measurements, and that the 
associated systematic errors are larger than the statistical errors recently discussed in the literature. 



I. INTRODUCTION 

A wide array of experimental strategies have been de- 
vised in order to identify the nature of dark matter (DM) 
P~H3] . A key parameter in many of these searches is the 
local density of DM, namely the density of DM particles 
in the solar neighbourhood, po- For instance, the rate of 
events in direct detection experiments, that seek to mea- 
sure the recoil energy in scattering events of DM particles 
off nuclei in the detector, is obviously proportional to the 
flux of DM particles through the detector, which in turn 
is directly proportional to the local DM density. Simi- 
larly, the neutrino flux from DM annihilation in the Sun 
is proportional to the capture rate of DM particles, in 
turn proportional to the flux of DM particles through 
the Sun, and therefore to p . As for indirect searches, 
the predicted flux of secondary particles, produced by 
the annihilation of DM particles, is proportional to pg. 
A careful determination of this quantity is therefore of 
paramount importance in order to extract the proper- 
ties of DM particles, especially when trying to perform a 
combined analysis of direct detection and LHC data [5j. 

Interestingly, as also pointed out in Ref . [6J , po is often 
assumed to be equal to 0.3 GeV cm~ 3 , with an error of a 
factor of 2 . However, this value is often given without a 
reference, and when a reference is given, it can be traced 
back to papers which are a few decades old, e.g. [3 [8]. 
Notice that the link between observed galaxy masses and 
their host dark matter halo masses has been studied ex- 
tensively using satellite kinematics HH] , weak lensing 
[TT] as well as via abundance matching [12] . Due to the 
large spread in halo characteristics (mass and concentra- 
tion) at a given galactic baryon mass, these studies are 
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not suitable for inferring detailed properties of our own 
Milky Way (MW) galaxy, e.g. the local dark matter den- 
sity at the solar radius. While the virial mass of the Milky 
Way is not tightly constrained (see discussion in section 
II A ) , more can be said about the local dark matter den- 



sity via mass modeling using observational constraints of 
the galactic baryons and halo stars e.g. [13] . 

A number of papers have appeared recently on this 
subject, where the authors attack the problem of the 
determination of po in light of recent observational re- 
sults [61 [HI EH- In Ref. dH (see also [15]). for instance, 
the authors considered a large set of observational con- 
straints of the Milky Way galaxy, e.g. the local stel- 
lar surface density, the angular rest-frame velocity at 
the solar radius and inner dark matter halo mass es- 
timates from the velocity dispersion of halo stars. By 
adopting a Bayesian approach to mass modelling of the 
Milky Way components, a local dark matter density of 
0.385±0.027 GeV/cm 3 (assuming an Einasto profile) was 
found. The quoted la errors are smaller than in previous 
studies due to the large set of input constraints, as well 
as the tight range that exists on a few of them (e.g. the 
combination of Oort's constants A — B). This value was 
found to be quite insensitive to the assumed spherical 
DM density profile. As we shall demonstrate, the non- 
sphericity of the dark halo and its reaction to galaxy 
formation introduces larger systematic uncertainties. An 
alternative technique has been proposed in Ref. |6j , where 
a constraint on p is obtained based on local observables 
and with presumably no dependence on the mass model 
of our Galaxy. 

Here we estimate the systematic uncertainty on po, 
with specific emphasis on the impact of departures from 
spherical symmetry. This has been studied in Ref. |16j 
using pure dark matter numerical simulations. We focus 
instead on a high-resolution simulation of a Milky Way 
like galaxy [T7J — which reproduces the correct proper- 
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ties of our Galaxy — and consider its realizations with 
and without baryons. In particular, extracting from the 
simulation with baryons the orientation of the Galactic 
stellar disk with respect to the DM distribution, we study 
the DM density for an observer located at ^8 kpc from 
the Galactic center on the stellar disk, and show that 
it is systematically larger than the average density in a 
spherical shell of the same radius. The latter is the ob- 
servable that has usually been inferred from dynamical 
constraints such as the local circular velocity, terminal 
velocities and velocity dispersions of tracer star popula- 
tions. Notice that we are disregarding fine grained struc- 
tures (such as microhalos and streams) since their effect 
on DM searches is likely negligible as shown in [15] . 

The paper is organised as follows: we introduce in sec- 
tion [n] the numerical simulations of Milky Way sized ha- 
los with and without baryons. In Section |TTT| we discuss 
the shape and enclosed mass of the DM component, and 
the impact on the determination of the local DM den- 
sity from dynamical measurements. We then conclude in 
Section HVl 



II. COSMOLOGICAL SIMULATIONS OF 
GALAXY FORMATION 

A. Models for the Milky Way halo 

In order to design a numerical model for the Milky 
Way halo, we first need to estimate our target halo mass 
and merging history. We then can explore large scale N- 
body simulations and pick a good Milky Way candidate, 
based on our a priori knowledge of the Milky Way halo. 
In this respect, it is usually admitted that a good Milky 
Way candidate should have its total mass between, say, 
5 x 10 11 M and 2 x 10 12 M PUDS], where the mass 
of a halo is defined here as M2000 the mass enclosed in 
a sphere of radius R200C in which the average density is 
200 times the critical density at that redshift. Although 
this mass range is relatively small (a factor of 4), choos- 
ing a large halo mass, say above 10 12 M , rather than a 
low halo mass, say below 10 12 M , has profound conse- 
quences on the baryons and dark matter dynamics and 
their final properties. This is due to the role played by 
the baryons inside the dark matter halo, and the fact 
that baryons dissipate their thermal energy, leading to 
the formation of dense and concentrated objects, usually 
in the form of centrifugally supported disks [20]. More- 
over, baryons are much better constrained by observa- 
tions than dark matter. For the Milky Way, we know 
that the disk can be decomposed in a pseudo-bulge of 
stellar mass Mb — 1 — 2 x 10 10 M and a rather ex- 
tended disk of mass M D -5-6x 10 10 M [HEE]. We 
also know that the Milky Way rotation curve is rather 
flat, with a maximum measured velocity of 220-250 km/s 
( [22] and references therein). For large halo mass models, 
this requires that the baryon fraction within the halo is a 
factor of 2 below the universal fraction (f2&/f2 m ~ 0.17), 



leading to the so-called "missing baryons problem" and 
requiring considerable amount of feedback to eject these 
baryons out of the halo boundaries ([33] and references 
therein). The second important consequence is that in 
these models, the halo is dark matter dominated. The 
dark matter distribution will be therefore only mildly 
affected by the baryons dynamics and pure N-body sim- 
ulations can be reliably used to estimate the local DM 
density at the Sun radius [21] ■ For small halo mass mod- 
els, on the other hand, the baryon fraction is close to the 
universal value. In this case, the effects of baryons dy- 
namics are maximal, with strong adiabatic contraction 
of the dark halo leading to a change of the global halo 
shape and concentration, as well as the formation of dark 
disks [2"51 - f2"8] . In this second scenario, in order to real- 
istically estimate the dark matter distribution, we need 
to model the baryonic physics and the formation of the 
stellar disk. We now describe such a model in the next 
section. 



B. Baryons dominated Milky Way model 

Agertz, Teyssier & Moore 17J (from now on ATM10) 
performed a large suite of cosmological simulations aimed 
to study the assembly of Milky Way like galaxies. The 
simulations were carried out using the Adaptive Mesh 
Refinement (AMR) code RAMSES, which includes treat- 
ment of dark matter, gas and stars. The gas dynamics is 
calculated using a second-order unsplit Godunov method, 
while collisionless particles (including stars) are evolved 
using the Particle-Mesh technique. The modelling in- 
cludes realistic recipes for star formation [29 , supernova 
feedback (SNII and SNIa), stellar mass loss, gas cool- 
ing/heating and metal enrichment [30] . 

ATM 10 demonstrated that realistic disk galaxies form 
if the efficiency of star formation is low at high redshift, 
due to small scale processes. We analyze the dark matter 
component in two of their hydro + TV-body simulations: 
SR6-n01elML and SR6-n01e5ML. The former simulation 
adopts an inefficient star formation law (1% per gas free 
fall time) and the latter a more efficient one (5%). These 
assumptions result in a strong difference in galactic disk 
size and concentration, hence dark matter contraction, as 
the bulge-to-disk ratio increases from ~ 0.25 to 1.3. This 
represents a shift in galactic Hubble type; the disk in the 
efficient simulation resembles an SO/Sa galaxy while the 
inefficient one is closer to a Sb/Sbc. These simulations 
were shown to bracket the observed Kennicutt-Schmidt 
relation from the THINGS survey [31] , and are therefore 
representative of typical z = disk galaxies. In addition, 
we compare the results from these simulations to a pure 
dark matter realization of the same halo. 

Note that these numerical simulations have a spatial 
resolution of approximately 340 pc (for the galaxies un- 
der study), which is good enough to resolve the bary- 
onic physics taking place in the stellar disk. Such is not 
achieved by other recent simulations that present resolu- 
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FIG. 1: The shape parameters b/a, c/a and T for SR6- 
nOlelML (solid black lines) and for the corresponding pure 
dark matter realization (dashed red lines), both at z — 0. Up- 
per (Lower) thin lines show b/a (c/a), whereas thick curves 
represent the triaxiality parameter T. Also shown are the 
Sun galactocentric distance Ro ~ 8 kpc and the virial radius 

Rvir = -R2OOC 



tions of a few kpc [3J1 [351 [33J [33]. EH and [37] studied 
the effect of disk formation on the local dark matter den- 
sity, specifically the formation of a dark disk component. 
These studies adopted a similar spatial resolution as in 
this work (few 100 pc), but the galaxies featured very 
massive bulge components, making them compare unfa- 
vorably to late type galaxies like the Milky Way. Hence, 
the lack of a massive bulge in our simulations signals a 
more realistic galaxy assembly, and possibly also a more 
realistic baryonic effect on the dark matter halo struc- 
ture. 



III. SYSTEMATIC UNCERTAINTIES ON THE 
LOCAL DARK MATTER DENSITY 

Any departure from spherical symmetry and any mod- 
ification of the DM enclosed mass due to adiabatic con- 
traction is expected to lead to a systematic error in the 
determination of the local dark matter density at the so- 
lar position, as in general it will be different from the 
average density on a shell of same radius, which is the 
quantity inferred from most dynamical measurements. In 
this section we study these two effects and quantify the 
corresponding systematic errors on the determination of 
Po- 



A. Halo shape 

It is well-known (e.g. [32]) that the inclusion of baryons 
in numerical simulations washes-out the prolateness of 
dark matter halos found in DM-only simulations. In or- 
der to measure the shape of the dark halo in the simula- 



tions under scrutiny (described in section |n]), we follow 
[25] [35] and compute for a given set of N p dark matter 
particles the matrix 



Jij — 



(1) 



where i,j = 1,2,3 index the coordinates in the reference 
system. The eigenvectors of Jij are the major, interme- 
diate and minor axes j a , % and j c , and the eigenvalues 
J a > Jb > Jc give the axis ratios through b/a = \J Jt/J a 
and c/a = \/ J c / J a - The major (minor) axis j a (j c ) corre- 
sponds to the axis around which the angular momentum 
is minimal (maximal). The triaxiality parameter 



T 



l-b 2 /a 2 
l-c 2 /a 2 



distinguishes prolate (T > 0.5) from oblate (T < 0.5) 
shapes. For a given R, we start by considering the par- 
ticles in the sphere of radius R to compute the principle 
axes, b/a and c/a. We then repeat the procedure select- 
ing particles in the ellipsoid u 2 + )2 + , 7 n)2 < R 2 , 



(b/aY ^ (claY 

where u, v and w are the coordinates along the major, 
intermediate and minor axes, respectively. The compu- 
tation is iterated until both b/a and c/a have varied less 
than 0.5%. 

In figure [I] we show b/a, c/a (thin lines) and T (thick 
lines) computed as described in the previous paragraph. 
Solid and dashed lines correspond respectively to SR6- 
nOlelML and the dark matter only realization of the cor- 
responding halo. As expected, in the absence of baryons 
the dark halo is manifestly prolate i.e. elongated along 
the major axis j a , while the numerical simulation with 
baryons produces a more oblate shape i.e. flattened along 
the minor axis j c . To check the orientation of the dark 
halo with respect to the baryonic component we plot in 
figure [i] the angle ip between j c and the normal to the 
stellar disk ri sc [- In SR6-n01elML the dark and bary- 
onic components are fairly aligned for R < 20 kpc. Fur- 
thermore, as the dot-dashed blue line indicates, in SR6- 
n01e5ML both components are even more aligned. No- 
tice that above ~100 kpc the presence of substructures 
affects significantly the shape measurement, as clear from 
figures [l] and [2] 

Now, we are interested in evaluating how the deter- 
mination of po in studies such as [14] is affected by the 
DM halo shape, in particular using the latest numerical 
simulations with baryons. As pointed out in |14j . local 

observables constrain efficiently 
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where v is the 
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circular velocity in our Galaxy and Ro the Sun galacto- 
centric distance. Such quantity depends on the baryonic 
content of the Milky Way and on the mass distribution 
of dark matter through the equations 
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FIG. 2: The angle between the normal to the stellar disk n s d 
and the minor axis j c . The solid black (dot-dashed blue) line 
corresponds to SR6-n01elML (SR6-n01e5ML) at z=0. The 
dashed red line shows the angle between the minor axis in 
the pure dark matter simulation and the normal to the stellar 
disk in SR6-n01elML. Also shown are the Sun galactocentric 
distance Ro ~ 8 kpc and the virial radius R V i r = ifeooc 
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where Kb encodes the contribution of baryons (at Ro 
mainly dominated by the disk component). In general, 
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where po is the spherically averaged dark matter density 
at i?o ~ 8 kpc. We aim at comparing the mean spheri- 
cal value po with the local one po in the MW-like simu- 
lated galaxy SR6-n01elML — such comparison yields the 
systematic uncertainty in the determination of the local 
dark matter density presented in works where spherical 
halos were assumed. In order to study this effect we 
select dark matter particles that lie inside the spheri- 
cal shell 7.5 < i?/kpc < 8.5. Given three orthogonal 
planes (e.g. the ones defined by the principal axes) we 
consider the portions of the shell lying at distances from 
each plane smaller than Aw/2 = 0.5 kpc. This proce- 
dure defines three orthogonal ring-like structures each of 
which we divide in equal parts encompassing an angle 
Aip = 7r/4 and thus a volume V = 2ir kpc 3 . For refer- 
ence in the following, 

10 7 M Q /kpc 3 = 0.38 GeV/cm 3 . 

Figure [3j sketches the dark matter density distribution 
in the above-mentioned rings for SR6-n01elML and the 
dark matter only simulation. In the former case we have 
used the stellar disk plane and two perpendicular planes, 
while in the latter the planes defined by the principle axes 
(of the set of particles in the shell 7.5 < i?/kpc < 8.5) 



were considered. In each plot the angle bins represent 
portions of the ring encompassing 7r/4 rad (as described 
before) and the vertical error bars are Poissonian. The 
sinusoidal-like modulations seen for both simulations are 
naturally expected due to the triaxiality of the halos. 
As a guiding line, we present in each plot of figure [3] the 
best fit function C\ +c 2 sin (2 (tp + C3)). Furthermore, one 
can appreciate large differences between the spherically 
averaged density po (dashed lines in figure [3]) and the 
density along each ring. 

We pick two extreme cases to bracket the systematic 
uncertainties in the determination of po- (i) the stellar 
disk plane in SR6-n01elML (figure [3| upper left panel), 
and (ii) the planes perpendicular to the minor and ma- 
jor axes in the dark matter only case (figure [3j lower left 
and lower central panels). Notice that we disregard the 
plane defined by the intermediate axis in the simulation 
without baryons because in that plane a stable baryonic 
disk cannot be formed 40J. In case (i), since the dark 
matter halo is flattened along the stellar disk, the local 
dark matter density is higher than the spherically aver- 
aged value: 

po/po = 1.01 - 1.41 . 

In case (ii) a broader range is obtained: 

po/po = 0.39 - 1.94 , 

in rough agreement with |16j . These values translate 
into systematic shifts on the local densities found e.g. in 
Ref. QJ]. 

For the sake of completeness a similar analysis was 
carried out for the simulated galaxy with efficient star 
formation rate, SR6-n01e5ML. For the equivalent to case 
(i) explained in the last paragraph, we obtain po/pa = 
1.21 — 1.60; note nevertheless that this extreme object is 
not MW-like. 



B. Enclosed mass 

Both theoretical arguments and simulations of galaxy- 
sized objects [351 133 [37] seem to indicate that baryons 
induce the contraction of the dark matter component to- 
wards the central part of the halo. This is also the case 
for the simulated galaxy SR6-n01elML when compared 
to the corresponding pure dark matter realization. Note 
that the dark matter only simulation needs to be rescaled 
down by a factor (f2 m — f^) /^m — 0.8333 in order to ac- 
count for the presence of baryons, and that one can only 
firmly trust the numerical simulation results above about 
2 times the resolution scale, i.e. 680 pc. We leave for a 
forthcoming work the detailed study of the DM profile 
and adiabatic contraction models. This topic, addressed 
in Refs. [Ml 123 [33] in the framework of different 
numerical simulations, is of particular importance for in- 
direct dark matter searches. 

A relevant quantity to analyse when trying to de- 
termine the local dark matter density is the enclosed 
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FIG. 3: The dark matter density in the spherical shell 7.5 < R/kpc < 8.5 along the stellar disk plane and two perpendicular 
planes for SR6-n01elML (top), and along the planes perpendicular to the principle axes for the pure dark matter simulation 
(bottom). The solid horizontal line represents the mean of the points and the dashed line shows the value of the mean density in 
the whole shell, dubbed po. The sinusoidal curve shown in each plot is the best fit to the points in the form ci +C2 sin (2 ((p + C3)). 



dark mass, M^ m {< R), since this affects the rota- 
tion curve of the Galaxy. Using the spherically av- 
eraged density reported in figure [3] (dashed lines), 
we find that the same enclosed mass Md m {< 8 kpc) 
would result in SR6-n01elML and the pure dark mat- 
ter halo if the local DM densities are rescaled such that 
p (SR6-n01elML)/p (DM only) ~ 0.9. The lower den- 
sity in the presence of baryons is simply a reflex of a 
more concentrated profile. In any case, these estimates 
do not translate directly into systematic uncertainties in 
the determination of po since precise determinations of 
local observables — namely the Oort's constants A ± B, 
the Sun galactocentric distance Rq and the local visi- 
ble matter surface dentity £*, see [HJ [T3] — constrain 



efficiently both po (through 



d{v 2 R) 
dR 



Ro 



) and M dm (< R ) 



(through v(Rq)). Therefore, we conclude that consider- 
ing a contracted dark matter profile would not change 
significantly the determination of the local dark matter 
density from precise dynamical observables, but would 
eventually prefer smaller concentration parameters (or, 
equivalently, larger scale radii R s ). 



IV. CONCLUSIONS 

The dark matter density in our neighbourhood is the 
key astrophysical ingredient that fixes the flux of DM 
particles crossing the Earth and the Sun, thus governing 
the scattering off nuclei in underground detectors as well 
as the capture rate in the Sun. Experiments looking for 
DM-induced nuclei recoils or neutrino fluxes from the Sun 
are hence crucially dependent on the local dark matter 



density. In the present work we have tried to quantify 
the systematic uncertainties associated to this parame- 
ter, and that affect determinations based on dynamical 
observables of our Galaxy. Using a very recent successful 
attempt to simulate a spiral galaxy that resembles the 
Milky Way, the dark matter density at the solar circle 
was analysed in detail and compared to the pure dark 
matter case. 

One major consequence of the inclusion of baryons is 
a significant flattening of the dark halo in the direction 
of the normal to the stellar disk, leading to a DM over- 
density in the local disk of up to 41% with respect to the 
spherically averaged value. More specifically, we found 
that in the MW-like simulated galaxy the local dark mat- 
ter density is higher than the spherically averaged value: 
Po/po = 1-01 — 1.41. In the DM-only case a broader range 
is obtained: po/po — 0.39 — 1.94. 

Ideally, one should repeat the analysis in Ref. [13] . i- e - 
a Bayesian approach to mass modeling of the Milky Way 
components, in presence of a triaxial profile like the one 
discussed here. However, based on the considerations 
presented above, a better estimate of the local dark mat- 
ter density can be obtained by raising by 21% the mean 
value obtained in Ref. [T3J for the spherical case, keep- 
ing relative statistical errors fixed and adding systematic 
errors. In the case of an Einasto profile, this procedure 
suggests 

pa = 0.466 ± 0.033(stat) ± 0.077(syst) GeV/cm 3 . 

Notice that the mean 21% enhancement with respect to 
the spherical local DM density is obtained for a specific 
simulated galaxy resembling the MW. The actual en- 
hancement in our Galaxy may, of course, be different, but 
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the main points here are that (i) the presence of baryons 
leads quite generally to a DM overdensity about the lo- 
cal disk, and (ii) the systematic uncertainties affecting 
Po are significant and, in some cases, already larger than 
the statistical ones. 

The baryons are also responsible for a non-negligible 
contraction of the DM distribution towards the central 
part of the galaxy. Even though this may be very im- 
portant in searching for products of DM annihilations — 
such as positrons or antiprotons from the galactic halo, 
and 7-rays or neutrinos from the Galactic Centre — wc 
found that it has no significant effect in the determina- 
tion of the local dark matter density using dynamical 
observables. 

An estimate of systematic uncertainties affecting the 
local dark matter density is an important step in assess- 
ing realistically our present knowledge on this key param- 
eter. Such knowledge is in turn an input in interpreting 



direct detection results, combining multiple DM-induced 
signals, and extracting compatible DM properties. An- 
other relevant ingredient for DM scattering and capture 
is the phase space density (see e.g. [551 155] ). namely the 
distribution of velocities f{\v |) — such topic is out of the 
scope of this work, but it would be interesting to study it 
in the set of simulations analysed here. In case of a pos- 
itive signal in direct detection experiments, for instance, 
the identification of the compatible particle physics pa- 
rameter space and the discrimination between different 
particle physics frameworks depend crucially upon the 
state of our knowledge in key astrophysical parameters 
such as the dark matter density in our neighbourhood. 
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